kcdf = function (x) 
{
    x = sort(x)
    n = length(x)
    den = distribution.default2(x, bw = "nrd0", kernel = "gaussian")

    if (n < 1) 
        stop("'x' must have 1 or more non-missing values!")
        
    rval <- approxfun(den$x, den$y, 
        method = "linear", yleft = 0, yright = 1, f = 0, ties = "ordered")
    class(rval) <- c("ecdf", "stepfun", class(rval))
    attr(rval, "call") <- sys.call()
    rval
}

distribution.default2 = 
function (x, bw = "nrd0", adjust = 1, kernel = c("gaussian", 
    "epanechnikov", "rectangular", "triangular", "biweight", 
    "cosine", "optcosine"), weights = NULL, window = kernel, 
    width, give.Rkern = FALSE, n = 512, from, to, cut = 3, na.rm = FALSE, 
    ...) 
{
    if (length(list(...))) 
        warning("non-matched further arguments are disregarded")
    if (!missing(window) && missing(kernel)) 
        kernel <- window
    kernel <- match.arg(kernel)
    if (give.Rkern) 
        return(switch(kernel, gaussian = 1/(2 * sqrt(pi)), rectangular = sqrt(3)/6, 
            triangular = sqrt(6)/9, epanechnikov = 3/(5 * sqrt(5)), 
            biweight = 5 * sqrt(7)/49, cosine = 3/4 * sqrt(1/3 - 
                2/pi^2), optcosine = sqrt(1 - 8/pi^2) * pi^2/16))
    if (!is.numeric(x)) 
        stop("argument 'x' must be numeric")
    name <- deparse(substitute(x))
    x <- as.vector(x)
    x.na <- is.na(x)
    if (any(x.na)) {
        if (na.rm) 
            x <- x[!x.na]
        else stop("'x' contains missing values")
    }
    N <- nx <- length(x)
    x.finite <- is.finite(x)
    if (any(!x.finite)) {
        x <- x[x.finite]
        nx <- length(x)
    }
    if (is.null(weights)) {
        weights <- rep.int(1/nx, nx)
        totMass <- nx/N
    }
    else {
        if (length(weights) != N) 
            stop("'x' and 'weights' have unequal length")
        if (!all(is.finite(weights))) 
            stop("'weights' must all be finite")
        if (any(weights < 0)) 
            stop("'weights' must not be negative")
        wsum <- sum(weights)
        if (any(!x.finite)) {
            weights <- weights[x.finite]
            totMass <- sum(weights)/wsum
        }
        else totMass <- 1
        if (!isTRUE(all.equal(1, wsum))) 
            warning("sum(weights) != 1  -- will not get true density")
    }
    n.user <- n
    n <- max(n, 512)
    if (n > 512) 
        n <- 2^ceiling(log2(n))
    if (missing(bw) && !missing(width)) {
        if (is.numeric(width)) {
            fac <- switch(kernel, gaussian = 4, rectangular = 2 * 
                sqrt(3), triangular = 2 * sqrt(6), epanechnikov = 2 * 
                sqrt(5), biweight = 2 * sqrt(7), cosine = 2/sqrt(1/3 - 
                2/pi^2), optcosine = 2/sqrt(1 - 8/pi^2))
            bw <- width/fac
        }
        if (is.character(width)) 
            bw <- width
    }
    if (is.character(bw)) {
        if (nx < 2) 
            stop("need at least 2 points to select a bandwidth automatically")
        bw <- switch(tolower(bw), nrd0 = bw.nrd0(x), nrd = bw.nrd(x), 
            ucv = bw.ucv(x), bcv = bw.bcv(x), sj = , `sj-ste` = bw.SJ(x, 
                method = "ste"), `sj-dpi` = bw.SJ(x, method = "dpi"), 
            stop("unknown bandwidth rule"))
    }
    if (!is.finite(bw)) 
        stop("non-finite 'bw'")
    bw <- adjust * bw
    if (bw <= 0) 
        stop("'bw' is not positive.")
    if (missing(from)) 
        from <- min(x) - cut * bw
    if (missing(to)) 
        to <- max(x) + cut * bw
    if (!is.finite(from)) 
        stop("non-finite 'from'")
    if (!is.finite(to)) 
        stop("non-finite 'to'")
    lo <- from - 4 * bw
    up <- to + 4 * bw
    y <- .C("massdist", x = as.double(x), xmass = as.double(weights), 
        nx = nx, xlo = as.double(lo), xhi = as.double(up), y = double(2 * 
            n), ny = as.integer(n), PACKAGE = "stats")$y * totMass
    kords <- seq.int(0, 2 * (up - lo), length.out = 2L * n)
    kords[(n + 2):(2 * n)] <- -kords[n:2]
    kords <- switch(kernel, gaussian = pnorm(-kords, sd = bw), # pnorm <- dnorm, -kords <- kords 
        rectangular = {
            a <- bw * sqrt(3)
            ifelse(abs(kords) < a, 0.5/a, 0)
        }, triangular = {
            a <- bw * sqrt(6)
            ax <- abs(kords)
            ifelse(ax < a, (1 - ax/a)/a, 0)
        }, epanechnikov = {
            a <- bw * sqrt(5)
            ax <- abs(kords)
            ifelse(ax < a, 3/4 * (1 - (ax/a)^2)/a, 0)
        }, biweight = {
            a <- bw * sqrt(7)
            ax <- abs(kords)
            ifelse(ax < a, 15/16 * (1 - (ax/a)^2)^2/a, 0)
        }, cosine = {
            a <- bw/sqrt(1/3 - 2/pi^2)
            ifelse(abs(kords) < a, (1 + cos(pi * kords/a))/(2 * 
                a), 0)
        }, optcosine = {
            a <- bw/sqrt(1 - 8/pi^2)
            ifelse(abs(kords) < a, pi/4 * cos(pi * kords/(2 * 
                a))/a, 0)
        })
    kords <- fft(fft(y) * Conj(fft(kords)), inverse = TRUE)
    kords <- pmax.int(0, Re(kords)[1L:n]/length(y))
    xords <- seq.int(lo, up, length.out = n)
    x <- seq.int(from, to, length.out = n.user)
    structure(list(x = x, y = approx(xords, kords, x)$y, bw = bw, 
        n = N, call = match.call(), data.name = name, has.na = FALSE), 
        class = "density")
}